!-------------------------------------------------------------LICENSE--------------------------------------------------------------!
!                                                                                                                                  !
!The MAP code is written in Fortran language for magnetohydrodynamics (MHD) calculation with the adaptive mesh refinement (AMR)    !
!and Message Passing Interface (MPI) parallelization.                                                                              !
!                                                                                                                                  !
!Copyright (C) 2012                                                                                                                !
!Ronglin Jiang                                                                                                                     !
!rljiang@ssc.net.cn                                                                                                                !
!585 Guoshoujing Road. Pudong, Shanghai, P.R.C. 201203                                                                             !
!                                                                                                                                  !
!This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License         !
!as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.             !
!                                                                                                                                  !
!This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of    !
!MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.                        !
!                                                                                                                                  !
!You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software     !
!Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.                                                   !
!                                                                                                                                  !
!-------------------------------------------------------------LICENSE--------------------------------------------------------------!

!==================================================================================================================================|
subroutine model_variables (ro, mx, my, mz, bx, by, bz, en, gx, gy, po, x, y, nx, ny)
!==================================================================================================================================|

   use parameters
   implicit none

   integer(4), intent(in) :: nx, ny

   real(8), dimension(nx), intent(in) :: x
   real(8), dimension(ny), intent(in) :: y
   real(8), dimension(nx, ny), intent(inout) :: ro, mx, my, mz, bx, by, bz, en, gx, gy, po

   real(8), dimension(nx, ny) :: te, pr
   real(8), dimension(0 : nx + 1) :: x_pre
   real(8), dimension(0 : ny + 1) :: y_pre
   real(8), dimension(0 : nx + 1, 0 : ny + 1) :: az_pre
   real(8), dimension(ny_total_finest) :: te_total_y

   integer(4) :: i, j, i_current, j_total, zeroy, e_dir
   integer(4) :: n_current

   real(8) :: amp_b, amp_v, d_current, l_current, start_x, start_y, ri2, rx_s, rx_e, theta_s, theta_e, r_s, r_e
   real(8) :: te_chr, te_cor, h_cor, w_tr
   real(8) :: g0, dx, dy, e_len, b_e, b_v, e_center, t_emerge_end, ri, r_x, r_y, r

!----------------------------------------------------------------------------------------------------------------------------------|
!  canopy
!----------------------------------------------------------------------------------------------------------------------------------|
   amp_b = 1.0d0

   n_current = 60
   d_current = 100.0d0
   l_current = 35.0d0

   start_x = -(n_current - 1) * d_current / 2 - l_current / 2
   start_y = -3.6

   g0 = -1.0
   te_chr = 0.6d0
   te_cor = 100.0d0
   h_cor = 8.5d0
   w_tr = 1.0d0

   if (.not. allocated(y_total)) allocate(y_total(ny_total_finest))
   if (.not. allocated(ro_total_y)) allocate(ro_total_y(ny_total_finest))

   if (ro_total_y_flag) then

      zeroy = (nint ((ny_total - 2 * ng) * zero_posy / rangey) + ng) * 2 ** (n_levels - 1) + 1
      dy = rangey / dble(ny_total - 2 * ng) / 2.0d0 ** (n_levels - 1)

      y_total(zeroy) = dy * 0.5d0

      do j = zeroy + 1, ny_total_finest
         y_total(j) = y_total(j - 1) + dy
      enddo
      do j = zeroy - 1, 1, -1
         y_total(j) = y_total(j + 1) - dy
      enddo

      do j = 1, ny_total_finest
         te_total_y(j) = te_chr + (te_cor - te_chr) / 2.0d0 * (tanh ((y_total(j) - h_cor) / w_tr) + 1.0d0 + eps)
      enddo

      ro_total_y(zeroy) = (2.0d0 * (te_total_y(zeroy - 1) + te_total_y(zeroy)) / 2.0d0 + dy * 0.5d0 * g0) /                        &
         (2.0d0 * te_total_y(zeroy) - dy * 0.5d0 * g0) * 1.0d3
      do j = zeroy + 1, ny_total_finest
         ro_total_y(j) = (2.0d0 * te_total_y(j - 1) + (y_total(j) - y_total(j - 1)) * g0) /                                        &
            (2.0d0 * te_total_y(j) - (y_total(j) - y_total(j - 1)) * g0) * ro_total_y(j - 1)
      enddo
      do j = zeroy - 1, 1, -1
         ro_total_y(j) = (2.0d0 * te_total_y(j + 1) - (y_total(j + 1) - y_total(j)) * g0) /                                        &
            (2.0d0 * te_total_y(j) + (y_total(j + 1) - y_total(j)) * g0) * ro_total_y(j + 1)
      enddo

      ro_total_y_flag = .false.

   endif

!   dx = x(2) - x(1)
!   dy = y(2) - y(1)
!   x_pre(1 : nx) = x
!   y_pre(1 : ny) = y
!   x_pre(0) = x_pre(1) - dx
!   x_pre(nx + 1) = x_pre(nx) + dx
!   y_pre(0) = y_pre(1) - dy
!   y_pre(ny + 1) = y_pre(ny) + dy

!   t_emerge_end = 30.0d0
!   e_center = 0.0d0
!   e_len = 5.0d0
!   b_e = 20.0d0
!   amp_v = 0.0d0

!   do j = 0, ny + 1
!   do i = 0, nx + 1
!      az_pre(i, j) = 0.0d0
!      if (y_pre(j) .le. 0.0d0 .and. abs (x_pre(i) - e_center) .le. e_len) then
!         if (t .le. t_emerge_end) then
!            az_pre(i, j) = b_e * cos (pi * (x_pre(i) - e_center) / e_len / 2.0d0) * t / t_emerge_end
!!         elseif (t .gt. t_emerge_end .and. t .le. t_emerge_end * 2) then
!!            az_pre(i, j) = b_e * cos (pi * (x_pre(i) - e_center) / e_len / 2.0d0) * (3 * t_emerge_end - t) / t_emerge_end / 2.d0
!         else
!!            az_pre(i, j) = b_e * cos (pi * (x_pre(i) - e_center) / e_len / 2.0d0) / 2.0d0
!            az_pre(i, j) = b_e * cos (pi * (x_pre(i) - e_center) / e_len / 2.0d0)
!         endif
!      endif
!   enddo
!   enddo

   t_emerge_end = 20.0d0
   e_dir = 1   ! 1: right to left; 2: left to right

   b_v = 0.0d0
   r = 8.0d0
   r_x = 0.0d0
   if (t .le. t_emerge_end) then
      r_y = -3.6d0 !- r * (1.0d0 - t / t_emerge_end)
      b_e = 36d0 * t / t_emerge_end
      !b_v = r / t_emerge_end * t / t_emerge_end * 5.0d0 !* 0.0d0
   else
      r_y = -3.6d0
      b_e = 36d0 * max(1.0d0 - t / t_emerge_end / 5.0d10, 0.0d0)
      !b_v = r / t_emerge_end * 5.0d0
   endif

   do j = 1, ny
   do i = 1, nx

      gx(i, j) = 0.0d0
      gy(i, j) = g0

      bx(i, j) = 0.0d0
      by(i, j) = 0.0d0
      bz(i, j) = 0.0d0

      te(i, j) = te_chr + (te_cor - te_chr) / 2.0d0 * (tanh((y(j) - h_cor) / w_tr) + 1.0d0 + eps)

      do j_total = 1, ny_total_finest
         if (y_total(j_total) .ge. y(j)) exit
      enddo

      if (j_total .eq. 1) then
         ro(i, j) = ro_total_y(j_total)
      else
         ro(i, j) = (ro_total_y(j_total - 1) * (y_total(j_total) - y(j)) + ro_total_y(j_total) * (y(j) - y_total(j_total - 1))) /  &
            (y_total(j_total) - y_total(j_total - 1))
      endif

      pr(i, j) = ro(i, j) * te(i, j)
      mx(i, j) = 0.0d0
      my(i, j) = 0.0d0

      ri = sqrt((x(i) - r_x) ** 2 + (y(j) - r_y) ** 2)
      if (ri .le. r) then
         my(i, j) = b_v * ro(i, j) * (r - ri) / r
         bx(i, j) = -b_e * e_dir * (y(j) - r_y) / (ri ** 1 + eps) !* (r - ri)
         by(i, j) =  b_e * e_dir * (x(i) - r_x) / (ri ** 1 + eps) !* (r - ri)
!         bx(i, j) = -b_e * e_dir * (y(j) - r_y) / (ri + eps)
!         by(i, j) =  b_e * e_dir * (x(i) - r_x) / (ri + eps)
      endif
      mz(i, j) = 0.0d0
   enddo
   enddo

   do i_current = 1, n_current

      rx_s = start_x + (i_current - 1) * d_current
      rx_e = rx_s + l_current

      do j = 1, ny
      do i = 1, nx
         r_s = sqrt((x(i) - rx_s) ** 2 + (y(j) - start_y) ** 2)
         r_e = sqrt((x(i) - rx_e) ** 2 + (y(j) - start_y) ** 2)
         theta_s = asin((y(j) - start_y) / r_s)
         theta_e = asin((y(j) - start_y) / r_e)

         if (x(i) .lt. rx_s) then
          theta_s = pi - theta_s
         endif
         if (x(i) .lt. rx_e) then
          theta_e = pi - theta_e
         endif

         bx(i, j) = bx(i, j) + amp_b * log(sin(theta_e) / sin(theta_s))
         by(i, j) = by(i, j) + amp_b * (theta_e - theta_s)
      enddo
      enddo

   enddo

   do j = 1, ny
   do i = 1, nx

!      bx(i, j) = bx(i, j) + (az_pre(i, j + 1) - az_pre(i, j - 1)) / 2.0d0 / dy
!      by(i, j) = by(i, j) + (az_pre(i - 1, j) - az_pre(i + 1, j)) / 2.0d0 / dx

      en(i, j) = (mx(i, j) * mx(i, j) + my(i, j) * my(i, j) + mz(i, j) * mz(i, j)) / ro(i, j) / 2.0d0 +                            &
         (bx(i, j) * bx(i, j) + by(i, j) * by(i, j) + bz(i, j) * bz(i, j)) / 2.0d0 + pr(i, j) / gmm1

      po(i, j) = 0.0d0

   enddo
   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
end subroutine model_variables
